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The nuclear-matter liquid-gas phase transition induces instabilities against finite-size density fluc- 
tuations. This has implications for both heavy-ion-collision and compact-star physics. In this paper, 
we study the clusterization properties of nuclear matter in a scenario of spinodal decomposition, 
comparing three different approaches: the quantal RPA, its semi-classical limit (Vlasov method), 
and a hydrodynamical framework. The predictions related to clusterization are qualitatively in good 
agreement varying the approach and the nuclear interaction. Nevertheless, it is shown that i) the 
quantum effects reduce the instability zone, and disfavor short-wavelength fluctuations; ii) large 
differences appear between the two semi-classical approaches, which correspond respectively to a 
collisionless (Vlasov) and local equilibrium description (hydrodynamics); iii) the isospin-distillation 
effect is stronger in the local equilibrium framework; iv) important variations between the predicted 
time-scales of cluster formation appear near the borders of the instability region. 

PACS numbers: 

Keywords: asymmetric nuclear matter, compact stars, liquid-gas spinodal instabilities, Skyrme forces, 
semi-classical theory, random-phase approximation 



I. INTRODUCTION 

In infinite nuclear matter, a phase transition of the liquid-gas type [H,IH is predicted to occur at sub-saturation den- 
sity with a critical temperature of about 15 MeV, depending on the theoretical models. This phenomenon corresponds 
to the separation of the system into two macroscopic (infinite) phases of different densities. This thermodynamic prop- 
erty is at the origin of finite-size spinodal instabilities, leading to the clusterization of the homogeneous medium. In 
this process, in contrast with the thermodynamic case, the Coulomb interaction can be included and density-gradient 
terms come into play: the liquid-gas instabilities can then be related to the formation of fragments occurring in finite 
nuclear systems. In heavy-ion collisions performed in nuclear facilities, around the Fermi energy, several signals of 
the liquid-gas phase transition have been analyzed @. In particular, the spinodal instabilities could be at the origin 
of the multifragmentation process [H [H, [|| . Several models have been built in order to describe the dynamics leading 
to the formation of clusters @, Q • Clustering is also expected to play an important role in compact stars 0, [l(| ■ It is 
predicted to lead to the presence of complex structures like the " pasta phases" occurring in the inner crust of neutron 
stars [O, [3 E3] • At lower densities and finite temperature, star matter may clusterize into light nucl ei fl5l [lt| . 
Around the core of supernovae, density fluctuations could have important effects on neutrino propagation [Tj], which 
is crucial for the description of supernova dynamics [l8| . 

This paper addresses the clusterization properties of an infinite nuclear medium: the onset of the liquid-gas instabil- 
ities is characterized by studying the incidence of a small density fluctuation introduced in the homogeneous medium. 
Such approach is based on a spinodal-decomposition scenario, which assumes that fragments are formed in an out- 
of-cquilibrium process. Indeed, if a homogeneous system is quenched into the region of finite-size instabilities, it will 
decompose into clusters: in the spinodal-decomposition scenario, the early dynamics governing this evolution consists 
in the rapid amplification of the unstable fluctuations, whose properties can be understood in the small-amplitude 
limit. For a fast reaction, fragments could be separated before equilibrium is realized, thus keeping the memory of the 
bulk instability. Identifying fragments produced through spinodal decomposition is then of great importance to obtain 
informations about the homogeneous-matter energetics at sub-saturation density. Experimental indications of this 
mechanism in heavy-ion collisions have been obtained for instance through the analysis of size correlations p^ . |2C| . 
We should note that the study of finite-size instabilities in homogeneous matter is complementary to the description 
of an equilibrated, clusterized system, as can be given at very low density by the virial equation of state [l5j]. With 
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TABLE I: Saturation properties of the three Skyrme forces SLy230a, SGII and RATP. Are given the values for the saturation 
density po: the binding energy E/A; the incompressibility modulus Koo; the symmetry energy a Sj o; the symmetry energy slope 
L = 3po(da B /dp); the symmetry energy curvature K sym = 9po(d 2 a B /dp 2 ); the critical temperature of the liquid-gas phase 
transition T c . 

this respect, the finite-size instability region we discuss can be viewed as the minimal region where the equilibrated 
system is formed of clusters. 

The model of infinite nuclear matter provides a useful insight into the behavior of nuclear systems, revealing robust 
qualitative features whose origin can be clearly analysed. The spinodal decomposition studied in this framework can 
be related to the nuclear multifragmentation process by neglecting the finite size of the homogeneous medium formed 
in the early stage of the collision. It also gives a picture of the instability against clusterization in compact-star 
matter, where the degree of freedom associated with the electrons can be neglected [211 ]. 

Our calculations are based on Skyrme density functionals. Specific attention is paid to the isospin degree of 
freedom, relevant for the physics of compact stars and exotic nuclei. The study of the instabilities is treated in 
three different approaches. Starting from the thermodynamic definition of the liquid-gas phase transition, we first 
concentrate on the free-ener gy c urvature properties. This leads us to study the finite-size density fluctuations in a 
Thomas- Fermi approach (2ll |22|. defining a static criterion for the existence of instabilities based on the free-energy 
variations. Afterward, consistently with this local equilibrium framework, we treat the time evolution by introducing 
hydrodynamical equations. This defines a first semi-classical approach of the clusterization dynamics. The second 
approach is fully quantal and dynamical: it is based on the linearization of the Time-Dependent Hartree-Fock theory 
(TDHF) [23|. This corresponds to the random-phase approximation (RPA). In this case, the evolution of the instability 
is described in the collisionless limit. The differences expected between hydrodynamics and RPA results have then 
two origins: the quantal effects and the collisionless description. In order to disentangle these two effects, we will 
also consider a third approach, the Vlasov method which is the semi-classical limit of the RPA. It then includes the 
collisionless description but not the quantal effects. 

Different Skyrme forces have been used in order to check the dependence of the results on the choice of a specific 
set of parameters. We take as a reference the Skyrme-Lyon parameterization SLy230a [24[, a modern force designed 
to treat neutron-rich matter, already employed for astrophysical applications (25|, [26[ . Two other Skyrme forces are 
used for comparison: RATP (27j - the first interaction which includes neutron matter into the fitting procedure - 
and SGII [28], a typical widely used parameterization. The interaction RATP is constrained by the neutron-matter 
equation of state calculated by Friedman and Pandharipande in the 80ies [29|, while the interaction SLy230a uses 
a more recent neutron-matter equation of state based on an improved realistic interaction [3(| • Let us remark that 
these interactions have quadratic velocity-dependent terms which modify the single-particle energies and give rise 
to effective masses m* (i=n, p). The effects of an isovector momentum-dependent mean field on nuclear-matter 
properties have been studied in Ref. [3l| . 

For these three interactions, Fig. [T] shows as a function of the density the binding energy E/A in symmetric nuclear 
matter and in pure neutron matter, as well as the symmetry energy a s . The different saturation properties are reported 
in Tab. |U One can see that the predictions of scalar properties, saturation density, energy and incompressibility are 
very close for these realistic forces. The values of the symmetry energy at saturation are also quite close, the main 
difference being its density dependence, as shown on Fig. [TJ The symmetry energy slope L and curvature K sym are 
also indicated in Tab. HI It should be noticed that the selected Skyrme forces present an "asy-soft" behaviour, namely 
values of L lower than 60 MeV, while recent constraints extracted from different experimental analysis indicate that L 
shall be larger than about 60 MeV [H, [H, [HI . This however has no drastic effects on the results we present, since the 
spinodal contour in asymmetric nuclear matter is non linearly related to the symmetry energy and its first and second 
derivatives, cf Eq. ([5]). A systematic analysis of the spinodal properties for various relativistic and non-relativistic 
interactions, "asy-soft" and "asy-stiff", is in preparation [35[. 

The paper is organized as follows. Before discussing the finite-size instabilities, we address the thermodynamic 
spinodal properties in Sect. [Til The hydrodynamical description for the clusterization process is treated in Sect. IIIII 
The collisionless dynamics based on the TDHF equation is treated in Sect. IIV1 where the quantal RPA and Vlasov 
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FIG. 1: (color online) Binding energy in pure neutron matter (top), in symmetric nuclear matter (middle) and symmetry 
energy (bottom) versus the density for the three Skyrme interactions considered in this paper (SLy230a, SGII and RATP). The 
saturation point in symmetric nuclear matter is also indicated (dots). 

methods are addressed. We discuss the different results in Sect.fV] and conclusions are drawn in Sect. I VII 

II. SPINODAL REGION FROM A THERMODYNAMIC AL APPROACH 

In this section, we address the thermodynamic liquid-gas phase transition of nuclear matter, where each phase is 
infinite and homogeneous. This allows a geometrical representation of the phase equilibrium, obtained by a Gibbs 
construction on the free-energy surface: thus can be emphasized the fundamental role of the free-energy curvature in 
the liquid-gas instabilities studied in the following sections. 

The thermodynamic equilibrium of a system is given by maximizing the entropy in the space of observables. 
For homogeneous nuclear matter, this surface presents an anomalous convexity in the phase-transition region. In the 
thermodynamic limit where interface contribution is neglected, this convexity can be corrected maximizing the entropy 
by linear interpolations which define pairs of phases in coexistence. In other words, the thermodynamic equilibrium 
is given by the Gibbs construction, which delimits the phase-transition boundaries. This coexistence region contains 
the spinodal region, where homogeneous matter is locally unstable (the entropy surface being locally convex). For a 
given value of the temperature, this is equivalent to a negative value of the free-energy curvature matrix, which we 
now determine. 

In the Skyrme density- functional theory (DFT), the average energy density of homogeneous nuclear matter TL = 
(H) /V is a function of the one-body density matrix. The full expression of Ti. obtained from a Skyrme-like interaction 
is given for example in Ref. 24]. Considering homogeneous, spin-saturated matter with no Coulomb interaction, the 
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energy-density functional reduces to: 



h 2 

n = 7i h = —t + C cSP t + D cSP3 t 3 
2m 

+C oP 2 + D pl + C 3 p a+2 + D 3P °p\ , (1) 

where the superscript "h" , introduced for later use, means "homogeneous" . Coefficients C an D are linear combinations 
of the standard Skyrme parameters [U,[33|- We have introduced the isoscalar and isovector particle densities, p and 
P3, as well as kinetic densities, r and t 3 : 

p = p n + P P ; r = r n + t p 

P3 = Pn - Pp i T 3 = T n - T p 

where, denoting i the third component of the isospin (n for neutrons and p for protons), the kinetic densities are 
defined by Ti = (p 2 )i/h 2 . The energy density ri is then a function of neutron and proton densities {pi,Ti}. From the 
Skyrme functional is derived the mean field for each type of particle i = (n, p). In uniform matter, the mean field 
reduces to a potential Ui and an effective mass m* such that the single-particle energies read: 

where p is the momentum quantum number of the free particles. In the grand-canonical formalism, individual levels are 
occupied according to the Fermi-Dirac statistics with occupation numbers ^(p) = [1 + exp(/3(ei(p) — , where 

[3 = l/T is the inverse temperature, and pn the chemical potential of the corresponding nucleon type i. The mean-field 
grand-canonical partition function Zq is formally identical to the free-Fermi-gas one, replacing pi by vi—pi — U, and 
m^ by m*. At given values of T, p n and p p , chemical potentials and kinetic densities are fixed, determining the energy 
density Ti. and the mean-field entropy density s, directly related to Zq. The free-energy density f(T, p n , p p ) is deduced 
from the Legendre transform i — H — Ts. We denote f h the free-energy density of a homogeneous system: 

t h (T,p n ,p p )=H h -Ts. (4) 

In the following, working at constant temperature, the dependence of the functions on T will be implicit for a more 
concise notation. For a given temperature, the thermodynamic spinodal region corresponds to a negative curvature 
of the surface i h (p n , Pp), «-e. a negative eigen-value of the curvature matrix C h = {C* }, with: 

whose lower eigen-value is denoted . 

Figure[2]shows the spinodal contours obtained for two values of temperature (T = and 10 MeV) and three Skyrme- 
like interactions: SLy230a, SGII and RATP. Small quantitative differences can be noticed between these interactions, 
especially at large asymmetry and density. These differences can be analyzed in terms of isovector properties of the 
density functional at (p n = p p ), which we will now discuss. We introduce the variables (p,y) where p is the total 
density and y = p 3 j ' p the asymmetry. Let us focus on the point (p = p s ,y = 0), defined as the high-density border of 
the spinodal for symmetric matter: the extension of the spinodal region in the asymmetry direction depends on the 
shape (convexity) of the spinodal contour at this point. This contour is defined by C\(p, y) = 0, where C< is the lower 
eigen-value of the free-energy curvature matrix. The isospin symmetry imposes [dC\/dy\p\ (p, 0) = 0. The convexity 
of the spinodal contour at the point (p s ,0) can then be characterized by the quantity C s = [d 2 C\/dy 2 \ p ] (p s , 0). 
Indeed, if C s is positive, for a small asymmetry y we have C< (p s ,y) > 0, meaning that the point (p s , y) is outside the 
spinodal contour. Thus, a positive (negative) value of C s corresponds to a concave (convexe) shape of the spinodal 
contour at point (p = p s ,y = 0). The value of C s depends on the symmetry energy and its density derivatives: 

C s = — [a s (p 2 a^' + 2pai) - 2{ P& ' s f] (6) 

where a^ = [da s /dp\ y ] (p = p s ,y = 0) and a" = [d 2 a, s /dp 2 \ y ] (p = p s ,y = 0). The respective contributions entering in 
Eq. ^ are shown in Tab. ILTlfor the three Skyrme-like interactions SLy230a, SGII and RATP. It is clear that none of 
these terms dominates the sum. Information on the density dependence of the symmetry energy (slope and curvature 
with respect to p) is then determinant to predict the qualitative behaviour of the spinodal. Conversely, information 
on the evolution of the spinodal extension with the asymmetry would bring constraints on the relation between these 
terms according to Eq. ©. 
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SLy230a 


0.1023 


25.73 


13.21 
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SGII 
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21.06 
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34.3 


RATP 


0.1026 


23.88 


12.69 


-12.41 


-20.2 



TABLE II: Contributions to the quantity C s coming from the symmetry energy and its derivatives at point p s (upper border of 
the spinodal in symmetric matter). The value of p s is also indicated. The Skyrme-like interactions SLy230a, SGII and RATP 
have been considered. 



III. LOCAL EQUILIBRIUM DESCRIPTION OF CLUSTER FORMATION 

Let us now turn to finite-size density fluctuations, relevant for the study of matter clusterization in neutron-star 
crust and supernova core, as well as nuclear multifragmentation. Coulomb and surface energy now contribute to reduce 
the instability induced by the thermodynamic spinodal. However, a region of finite-size instabilities can still be found, 
causing the instability of homogeneous matter against cluster formation. This region is determined considering plane- 
wave density fluctuations, varying the corresponding wave-number q: indeed, any density fluctuation can be expressed 
as a coherent sum of plane waves and the instability with respect to at least one q value signs the presence of finite-size 
instabilities. 

This section presents a first approach for the study of the unstable region: we assume local equilibrium in a semi- 
classical framework, the local densities defining local Fermi momenta. First, an instability criterion is defined from 
a Thomas-Fermi calculation of the free-energy variation; then hydrodynamical equations are introduced in order to 
describe the temporal evolution. 
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SLy230a 


84.5632 


70.8227 


77.693 


6.87028 


SGII 


43.275 


66.3604 


54.8177 


-11.5427 


RATP 


80.0409 


79.9703 


80.0056 


0.0353125 



TABLE III: Coefficients of the density-gradient terms in the Skyrme density functional. 



A. Finite-size instabilities in the Thomas-Fermi approach 



We detail here how the free-energy variation induced by a small density fluctuation is derived from the Thomas- 
Fermi approximation. Let us consider a plane- wave density fluctuation of wave number q, affecting both the neutron 
and proton densities (i = n, p): p, = po,i + $Pi, with Spi — A;e lqr + cc. The size of a density fluctuation, A, is 
characterized by the momentum q: X — 2n/q (note that the thermodynamic limit is obtained back with the infinite- 
wavelength density fluctuation, namely q — > 0). The Thomas- Fermi approximation supposes that density variations 
6 pi are smooth enough to allow the definition of a local Fermi momentum. This approximation is valid if the number 
of states piQ, - where fi is the volume - is very large for a typical scale of A. Taking for the volume £1 = A 3 , this 
condition leads to the following criterium: 47rg(fcF i /q) 3 /3 3> 1, where g is the spin degeneracy g = 2. It is then clear 
that this approach is valid if q is smaller than the Fermi momentum fcFi @> HH • 

In the description of nuclear matter with density fluctuations, finite-size effects appear because of the momentum 
dependence of the nuclear force, which is also related to its range. In the Skyrme framework, the surface energy can 
be written as: 

U v = C n v n (V Pn ) 2 + C p v p (V Pp ) 2 + 2C n v p V Pn • V Pp , (7) 
where coefficients Cy are combinations of the standard Skyrme parameters: 

Cj n = = | Ml - *i) - *a(l + x 2 )] , rsl 
Cn V P = Cj n = | [3t x (2 + Xl ) - t 2 (2 + x 2 )] . W 

To get a deeper insight, we can also express 7i v in terms of the isoscalar and isovector density gradients, namely: 

n v = c£(Vp) a + c£(Vp a ) 2 (9) 

where the isoscalar and isovector coefficients are related to (171 as: 



fiV f<V f~<V i r~<V 

u nn — ^pp — °11 T <^33 i 

CV Wv /~>V 
np ^ °pn ^ °11 "^SS 



(10) 



The values of the coefficients C v associated with the density-gradient dependence of the Skyrme functional are given 
in Tab. IIIII for several interactions. We can notice that the main part of the surface energy is due to the isoscalar 
density-gradient term, associated with the coefficient C^, which is much smaller for SGII than for SLy230a and 
RATP: consequently, large density gradients should cost less energy with SGII. 

The local- Fermi- momentum approximation generates a coordinate-dependent Fermi momentum kp(r). Together 
with the temperature T, it provides expressions for the local particle, kinetic-energy and entropy densities. Using the 
Skyrme energy functional and the entropy, we can directly evaluate the free-energy density, which can be decomposed 
into a bulk and a surface term. 

The bulk term contains the bulk part of the energy functional (Eq. ([T])) and the entropy term. It is obtained from 



the space average of the free energy density f h (defined in Eq. calculated with the local densities 2l|. Performing 
a development up to the second order in the amplitudes A$ associated with the density fluctuations - the first order 
being null due to particle conservation - the bulk variation of the free-energy is given by: 

c f b _ \ " A ' A i + ■fyido.n, A), P ) ni , 

^ 2 dp, ' { ' 

which depends on the curvature-matrix elements of the homogeneous free-energy C^- = dpi/dpj. 
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The surface term is generated by the gradient term of the energy functional. The space average of TL V leads to the 
surface contribution: 

S£V = TiJ ^ V(r)dr = q2 ? (A ^* + A * Aj) c ^ • (12) 

Considering finite-size fluctuations, one shall also take into account the effect of the Coulomb interaction. This was 
not the case in the thermodynamic framework discussed in Sec. HH since it is well known that the Coulomb energy 
diverges in infinite charged matter. In star matter, electrons equilibrate the positive charges and remove the divergence. 
In finite nuclei, the Coulomb energy might be large but finite. In both cases, the physically important quantity is 
not the total Coulomb energy but its variation S£ c induced by density fluctuation. In the present calculation, we 
have neglected the small contribution of the exchange Coulomb interaction so that the Coulomb contribution to the 
free-energy variation reads: 



A-KCi^ \ - A t A* + A*A 7 - 
q 2 



i.j 



where e p = g e /v47reo, e n = 0, with q Q the elementary electric charge. 

Finally, the total free-energy variation is given by the summation of the three terms, 5i = <5f b + <5£ v + 5£ c : 



s{ AjAj + A* A, ^ i{po Po , p ) + + ^ep\ 

XI 2 V <>/', <r J 

which can be expressed in a matricial form in the space of the density fluctuations A = (A n , A p ): Si ~ A*C f A, where 

c '^- ! (SS)v(s:)- 

The curvature matrix C f is a function of the temperature T, the average densities and the wave number of the 
fluctuation q. It contains the bulk term of the curvature, C h , the density-gradient term of the nuclear interaction 
proportional to q 2 , and the Coulomb- interaction term proportional to 1/q 2 . These two additional terms being positive 
tend to reduce the instability induced by the bulk term. Thus, the instability will be quenched at large values of q 
by the surface contribution, and at small values of q by the Coulomb interaction. It is then interesting to evaluate 
the curvature matrix for intermediate values of q, where it may still acquire negative values. In the following, the 
calculations show that there is still a region where the lower eigen-value of C f is negative for a finite interval of q 
values, leading to the development of spinodal instabilities. 



B. Hydrodynamical evolution 

The free-energy curvature analysis presented so far provides a static description of the finite-size instabilities. In 
order to describe the early dynamics of matter clusterization induced by such instabilities, we use a hydrodynamical 
framework, which corresponds to the local equilibrium limit and so is consistent with the above Thomas-Fermi 
approximation of (5f. Let us first consider a single fluid described by the Euler and continuity equations [7|: 

d(mpv)/dt = — VP , (16) 
dp/dt + V(H = , (17) 

where P is the system's pressure, m is the particle mass and p the particle density. In the limit of small-amplitude 
fluctuations p — po + Sp, Eqs. (|16tilT[) can be transformed into a diffusion equation: 

d 2 p 1 dP(pp) , po d 2 i(p ) 2 

~ a V p — —i^— V p . (1») 

or m dp m dp 

which is obtained using the relation between the pressure and the free-energy density f, dP/dp = pd 2 f/dp 2 . Expanding 
Eq. (|18[) on a plane-wave basis gives the following dispersion relation for the density modes: 

^ = Po d^po) q2 
m dp 
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If the free-energy curvature d 2 f/dp 2 is negative, the frequency is purely imaginary, u> = ry, and induces an exponential 
growth of the fluctuation with typical time r = 1/7. In asymmetric unstable nuclear matter, we generalize Eq. (|19[) 
introducing fluctuations along the unstable direction of the curvature matrix. This leads to : 

- 2 = -7 2 = ^q 2 , (20) 
m I 

where m is the nucleon mass, p the total average density and C< the lower eigen-value of the free-energy curvature 
matrix in the space of density fluctuations. The factor 1/2 is due to the change of variables from (p, P3) to (p n ,p p ). 
Equation (|20p characterizes the properties of the finite-size instabilities in terms of the free-energy curvature matrix, 
in the framework of a semi-classical and local equilibrium approach. In the following, it will be denoted as the 
hydrodynamical approach {4o| . 



IV. COLLISIONLESS DYNAMICS OF CLUSTER FORMATION 



Let us now consider a fully quantal approach for the description of the growth of density instabilities. This approach 
is based on the linear response theory to unstable fluctuations. In the hydrodynamical approach, a local equilibrium 
was assumed: in microscopic words, it means that for each particle species a Fermi sphere centered on zero momentum 
was defined in each point of space. Now, in contrast, each particle species is associated with a single Fermi sphere 
affected by a momentum shift (particle-hole excitations) causing density fluctuations in coordinate space. These two 
limits correspond respectively to first and zero sound. 

The linear response function is obtained from the small-amplitude limit of the TDHF theory (Time-Dependent 
Hartree-Fock): this corresponds to the RPA (Random-Phase Approximation) [23j]. In the present work we use the 
techniques developed in Refs. [H, 39| to obtain the RPA response at finite temperature, in the peculiar case of Skyrme- 
type effective interactions. We first solve the RPA in a fully quantal framework. We then introduce the semi-classical 
limit of the RPA, namely the Vlasov approach, in order to identify the quantum effects in the comparison between 
hydrodynamics and RPA results. 



A. Quantal effects described in the HF+RPA framework 

In the quantal framework, we have to linearize the TDHF equation, which describes the time evolution of the one- 
body density matrix pij = (ctcj), cj; being the annihilation (creation) operator associated with the single-particle 
state i: 

i^ = [h + V ext ,p] , (21) 

where h is the mean-field operator and V ex t is an external field. In the absence of perturbation, the mean-field ground 
state is given by po and h , which is the Hartree-Fock solution satisfying [h ,/3 ] = 0. In the following, we consider 
spin-saturated systems. For infinite nuclear matter, the operators ho and po are diagonal in momentum space and we 
have: 

f p Q |k> = sff k |k) 

\ h |k) = ^ k |k) 

where g is the spin degeneracy, and fk the Fermi-Dirac occupation for a level of momentum k. 

In this section, we consider a large system of volume ft and the plane waves |k) are defined such that (r|k) = 
e ik ' r /Vn. 



1. Linear response to an external probe 



We calculate the response function of nuclear matter at finite temperature to an infinitesimal external field of the 
form: 

V cxt - £ e^-iM-fcj)* (t < 0) ; V cxt =0 (t > 0) , (23) 
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where r\ is a positive infinitesimal part ensuring the adiabatic switching of the external field. The excitation carries 
a momentum q and an energy u>. For simplicity, we first treat nucleons as a single-component fluid. At late enough 
times and small enough external field £, the expectation value of the density fluctuation has the same space and time 
dependence as the external field: 

Sp(r, t) = aji-'-Ku+W (24) 
The linearization of Eq. (f2Tj) leads to the polarization function: 

tt / s ol n (u,q) 

n(w ' q) = ^ = i-n (., q )v res - ^ 

In this expression, we have introduced the Lindhard function Ilo (corresponding to the response function in the 
absence of residual interaction): 

n (<J,q) = g „ /dk— — ^±3_ ; w kq =w k+q -w k , (26) 

and the residual interaction v rcs which accounts for the variation of the mean field <5h due to the density fluctuation. 
Details for the derivation of the response function (|25p are given in Appendix [X] 

A collective mode (u>, q) of the nuclear fluid corresponds to a large enhancement of the response function, or a small 
value of the denominator of Eq. ([25]). In the spinodal unstable zone, the response function diverges for a given value 
of to and q. The dispersion relation w(q) for the unstable mode is then defined by (1 — IIo(a;(q), q)v rcs = 0). 

Let us now take into account the isospin degree of freedom. In this case, we distinguish neutron and proton fluids, 
and we get the response function and the residual interaction in terms of the matrices: 

n = ( n °- no ° ) , (27) 



v res v res 

nn np 

v res v res 

pn pp 



PP 



(28) 



n = (i - nov-T 1 n . (29) 

The linear response II(w, q) presents a divergence for the zero of matrix (1 — IIoV res ), corresponding to the condition: 

det (1 - n v rcs ) = . (30) 



2. Residual interaction and free-energy curvature 

We now have to express the residual interaction v res entering the definition of the response function (|129"|) . We give 
in Appendix [B] the steps demonstrating the link between v rcs and the free-energy curvature matrix C discussed in 
Sect. IIIII To summarize, considering a density fluctuation 

<5Pi( r ) = ^2 Phi ; /'-.q A,. q .' l<i r (i = n,p), (31) 
q 

we have the average residual-energy variation per unit volume: 

^ rcs = EE^^W (U = n,p)- (32) 
q i,j 

The following expression is obtained for the residual-interaction elements as a function of the free-energy curvature: 

<f(q) = C^(q)-^, (33) 
where No ; i is the density of states at Fermi level for the particle species i. This corresponds to the matrix expression: 

C f = (No)" 1 + ; No = ( N J n N ° ) . (34) 
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We remind that according to the semi-classical free-energy criterion discussed in the previous section, an instability 
occurs when the matrix C f has a negative eigen-value. The corresponding spinodal region is circled by a frontier for 
which the lower eigen-value is zero, responding to the condition: 



detC f = dct 



(No)" 







(35) 



Now, the RPA modes arc defined by the condition (|30[) . which can be expressed in a closely related form as: 



detC RPA = det 



(No)" 



Pov 1 



= 



where 



P (w,q)=-(N ) 'no^q) 



(36) 



(37) 



Note that the matrix C RPA which has just been introduced, in contrast with any curvature matrix, is not symmetric 
in the general case. 

In the RPA framework, the spinodal frontier is defined by the occurrence of a zero transfered energy mode. For 
T = 0, in the limit q — > the uncorrelated response function then reduces to [40l |: 



lim n . 4 (w = 0, q) = -N ,., 



(38) 



According to Eqs. (|35H37|) , this implies that, for T = and q — ► 0, the condition defining the RPA spinodal frontier 
reduces to the semi-classical one, det C f = 0. Note however that, taking into account the Coulomb interaction, 
the instability interval of q has a lower limit, so that the RPA and semi-classical spinodal frontiers always remain 
separated. 



3. Unstable modes in the RPA framework 



In the RPA approach, the instability region is characterized by the occurrence of imaginary modes: u) = w T + ry. 
The imaginary term i7 leads to an exponential increase of the fluctuation amplitude with the typical time r = I/7. 
Considering the general expression of u in the complex plane, the uncorrelated response function expressed by Eq. 
can be written (for each particle species i): 



no,t(w,q) 



h(2nf 



dk 



f; 



I. 



iM+q 



2<? 



dkf, 



(39) 



'kq 



The condition ([30]) can be satisfied only if the imaginary part of IIo cancels, which according to Eq. ([39| implies that 
lu 2 is real. Consequently, unstable modes correspond to a purely imaginary energy such that u = vy. The resulting 
response function reads: 



IIo,i(7,q) = -- 



2.9 



dkf 



i.k 



t^kq 



7 2 +^q 



At zero temperature, this integral has an analytic solution: 

n o,*(7,q) = -N , 4 Po, l (7,q) 

where, introducing the quantity wf = 7/[S/2m*] (wi has the dimension of a momentum), Po,i is given by: 



(40) 



(41) 



Po,j(7,q) = PoAmh),q) 



1 



2A:Fiq J 



kdk In 



(2fcq + q 2 ) 2 + w 



(2/fcq-q 2 ) 5 



Introducing the reduced variables i>i = wf /2qkp i = m*7/(qfcpi) and qF ; = q/2fcF i , we obtain the analytical expression: 



p o, 4 (7,q) = p o, 4 (^(7 I q) I qF,(q)) 



i + ^-q| 41 (i + q Fl ) 2 + ^ 2 



4qp 4 



■In 



(i-q Fl ) 2 + ^ 



arctan 



1 + qF, 



arctan 



1 - qFi 



(42) 
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which verifies, consistently with Eq. (|38)) . the condition at the spinodal frontier: 

limP 0>i (7 = 0,q) = l. (43) 

The RPA dispersion relation for the unstable mode 7(q) is deduced from Eqs. (|36|) and (|42|) . In this quantal 
approach, for given temperature and densities, uniform matter is unstable against a fluctuation of finite momentum 
q if detC^ A (7) presents a zero for a given value of the typical time r = 1/7. The corresponding vector gives the 
direction of the growing fluctuation in the density plane. 

At finite temperature, the same procedure can be followed but the response function defined by Eq. (I40| needs 
numerical calculation. Two different methods have been considered: one involves the calculation of Fermi integrals (4(il | 
and the other consists in performing a pondered sum of zero temperature response functions [381 ] : 

Re n 0ji (7, q, T) = - /Re n 0>i (7, q, T = 0, k F = k)df t (k, T) . (44) 

Jk 

At zero temperature, df^(fc) = — 5(k — k Fi )dk, yielding the well-known expression for RelLj^. Both methods have 
been used in order to check the accuracy of our results at finite temperature. 



B. Vlasov approach 

We now introduce the semi-classical limit of the RPA. The Vlasov approach describes the matter response to a 
momentum transfer q and energy transfer to like the RPA approach, but in a semi-classical approximation. It is 
then an intermediate description between the hydrodynamics and the full quantal RPA. Again, for simplicity, we give 
expressions for the case of a single fluid. The semi-classical version of the TDHF equation is obtained in the limit 
h — > [HI]. Using a Wigner transform, the commutator divided by ih is then replaced by a Poisson bracket, where 
the one-body density and potential operators are expressed by functions of time, position and momentum (t, r, k). 
This gives the Vlasov equation, describing the semi-classical time evolution in the collisionless limit: 



df(*,r,k) 
dt 



{h(t,r,k)+V ext (t,r,k) , f(i,r,k)} 



(45) 



where f is the occupation function in phase space, h the mean field and V cxt an external potential. In the unperturbed 
state, the mean field is h (k) and the occupation at temperature T — 1/(3 is given by the Fermi-Dirac distribution 
fo(k) = [1 + exp(/3(ho(k) — pi))]. The density fluctuation corresponds to the occupation variation i(t, r, k) = fo(k) + 
5i(t, r, k), inducing the mean-field variation h(t, r, k) = ho(k) + 5h(t, r, k). 

Let us stress here the difference between Thomas- Fermi and Vlasov approaches. Both approaches give a semi- 
classical description of the state occupations, using the function i(t, r, k): however, they define the perturbation 
<5f(i, r, k) in different ways. In the Thomas-Fermi approach, the variation of the local density is associated with a 
variation of the chemical potential, according to the local equilibrium description, and the occupation function is 
given by a Fermi-Dirac distribution: f(i,r,k) = [1 + exp(/3(ho(k) — /x(t, r)))]. In terms of the density matrix, this 
corresponds to a variation of the Fermi surface. In contrast, in the Vlasov formalism f(t, r, k) is given by a canonical 
transformation of fo(k): the variation of the distribution can be expressed in terms of a generating function S(t, r, k) 
such that Si = {S,fo} [H|]. In terms of the density matrix, this corresponds to the excitation of particle-hole states 
with a finite momentum transfer. 

The linearized Vlasov equation (semi-classical limit of the RPA) can be solved following a method similar to the 
one exposed in Ref. (44j . The Vlasov eigen-modes for asymmetric nuclear matter obey the matrix relation: 



where (with i = n, p): 



L(s) 



detC vlasov (w,q) =det 



(No)" 



Lv 1 



_ ( L(s n ) 
L(s F 



dh j({fc, = k Fl }) 
1 dk 



= 0, 



L(si) = 2 - Si In 



(46) 



(47) 



L(sj) denotes the semi-classical Lindhard function. The low-q limit of the uncorrelated response function no^c^q) 
is related to L(si) in such way that we get in terms of matrixes: 



lim(-N n (w,q)) = lim P (w,q) = L/2 

q^O q— >0 



(48) 
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FIG. 3: (color online) Dispersion relation at a given density and asymmetry, for zero temperature. Dots represent the boundaries 
of the unstable interval of q. Top: 1/r is given as a function of the wave number q. Bottom: associated directions of phase 
separation in the density plane Sps/Sp. The horizontal line at 0.4 is the direction of constant proton fraction pzj p. Left: 
three approaches applied with SLy230a: RPA (thick line), Vlasov (dashed line), hydrodynamics (thin line). Right: RPA and 
hydrodynamics with different Skyrme forces. 



This means that the RPA eigen-modes obeying Eq. (|36|) reduce to the Vlasov ones in the low-q limit (whatever 
the value of u>). On the other hand, for u = we have L(0)/2 = 1, so that C vlasov (w = 0,q) reduces to the free- 
energy curvature matrix C f (q): consequently, the frontier of the spinodal region is the same for the two semi-classical 
approaches. On this common frontier, the Vlasov and hydrodynamics eigen-modes are exactly identical. Differences 
shall appear inside the spinodal region where lo ^ 0. 



V. DISCUSSION OF THE RESULTS 



In this section, we present the properties of matter clusterization obtained from the different approaches presented 
above. The quantal RPA is compared to two semi-classical approaches: Vlasov and hydrodynamics, corresponding 
respectively to the collisionless (zero sound) and local equilibrium (first sound) limits. We first consider instability 
against density fluctuations for several momenta q, exploring the dependence on the temperature and on the nuclear 
interaction. Then, we concentrate on the most unstable mode defined as the mode for which the typical instability 
time is minimum and explore the clusterization properties. 



A. Unstable eigen-modes 



Let us determine the different instability properties associated with a given wave-number q. The unstable eigen- 
modes obey the dispersion relations illustrated for fixed densities on Fig. [3] The upper part of this figure represents the 
inverse typical time l/r(q), which characterizes the time evolution of the modes. The quantum effects can be directly 
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appreciated from the comparison between RPA and Vlasov results. Although these results are quasi-similar at the 
low-q limit, as expected, they separate for higher values of q. We then observe a quantum quenching of the instability, 
as discussed in Ref. [36] : the RPA unstable modes are associated with lower values of 1/r (i.e. a slower amplification of 
the fluctuation) and they end at a lower value of q. Indeed, in the Vlasov approach the uncorrelated response function 
is approximated by its limit at low q, which amounts to neglecting a part of the kinetic-energy cost of the fluctuation. 
In the RPA where the exact value of Hq(uj, q) is calculated, the high-q density fluctuations are disfavored. Let us now 
consider the semi-classical approach based on local equilibrium (hydrodynamics). It presents the same qualitative 
features as the previous ones: in particular, as discussed earlier, the instability region corresponding to a negative 
value of the free-energy curvature is formally the same as for the Vlasov approach. However, the hydrodynamical 
dispersion relation calculated from the curvature C< leads to a much faster evolution, 1/t being typically twice as 
large as for the RPA. It also presents a less symmetric behavior, with a sharper slope toward the high-q border: this 
is due to the relation 1/r 2 cx |C<|q 2 , in which the factor q 2 enhances the high-q side of the dispersion relation. 

The lower part of Fig. [3] deals with the direction of phase separation, also included in the dispersion relation 
determining the eigen-modes. This feature is connected with the isospin distillation, which is a well-known aspect 
of the liquid-gas phase transition in asymmetric nuclear matter [42j : in order to minimize the symmetry energy, the 
dense phase is more symmetric than the dilute one. This is observed in the construction of the thermodynamic phase 
equilibrium, and reflected by the spinodal instability direction: instead of following a line of constant composition in 
the density plane, phase separation is rotated towards the isoscalar direction. For a neutron rich point of densities 
(p, p 3 ), for which phase separation direction is given by the vector coordinates (5p, Sp 3 ), the normal isospin distillation 
is reflected in the relation: 

< 5p 3 /Sp < p 3 /p (49) 

We have checked that this relation is verified by the thermodynamic spinodal instability for all the Skyrme-like 
interactions we consider in this article. Moreover, it has been shown in Refs. [43| that for a set of realistic Skyrme- 
like and Gogny interactions, the instability is still isoscalar-like, which means that 5p 3 /8p < 0.5 in a domain of 
reasonable isospin asymmetries. Now when finite-size instabilities are considered and the electric charge comes into 
play, the Coulomb interaction disfavors proton-density fluctuations: this goes against the normal distillation effect 
(for neutron-rich matter). As can be seen on Fig. [3l isospin distillation is particularly affected at low values of q, 
where Coulomb effects dominate: an inverse distillation is even found near the low-q border of the instability interval. 
However, the normal isospin distillation remains the dominant feature, especially if we consider the most unstable 
modes [2ll . |4J] : it leads to the formation of clusters more symmetric than the homogeneous initial system, immersed 
in a more neutron-rich gas. 

Figure [3] shows a very close agreement between Vlasov and RPA results as for the value of Sp 3 /6p. The direction 
corresponding to the hydrodynamical approach is nothing but the direction of minimal free-energy curvature, i.e. 
the eigen- vector u< associated with the lower eigen- value of the curvature matrix C f . Qualitatively, it presents a 
similar evolution; furthermore, once again hydrodynamics and Vlasov curves are constrained to join at the borders 
of their common instability interval. However, inside this interval the hydrodynamics curve is clearly distinct from 
the two collisionless approaches. This shows that the phase separation direction followed in a collisionlcss approach 
is quantitatively different from the direction u< of the negative free-energy curvature: it leads to a slightly reduced 
distillation effect. 

The right column of Fig. [3] shows the dependence of the dispersion relation on the Skyrme force. RPA and 
hydrodynamics results are given for SLy230a, SGII and RATP. All present similar features, ensuring that the results 
elsewhere shown with SLy230a alone are qualitatively representative. The quantitative differences between SLy230a 
and RATP are weak. SGII however favors modes of larger q-values, and presents an extended instability interval on 
the high-q side. This difference comes from the surface properties of the forces: as can be seen in Tab. IIII1 the energy 
cost for the introduction of a density gradient is lower with SGII. 

Let us now address the extension of the instability in the density plane. We first consider instability regions 
associated with given values of the transfered- momentum q. In a condensed way, we name such regions " q-spinodals" . 
Inside a q-spinodal contour, a density fluctuation of momentum q is unstable and amplified according to the dispersion 
relation discussed above. Within the RPA approach, the q-spinodal contour at temperature T is determined by the 
divergence of the linear response H(lj = 0, q, T). In the semi-classical description (for both Vlasov and hydrodynamics), 
the q-spinodal contour corresponds to a zero eigen- value of C f (q, T). As shown by Eq. both curves are equivalent 
at the limit q — > 0. Figure 0] represents different RPA and semi-classical q-spinodals at T = 0. For q = 50 MeV/c 
almost no difference can be seen in the density plane. On the lower part of the figure, the q-spinodals are drawn 
in the plane of Fermi momenta, k-p i — (Sir 2 Pi) 1 ^ 3 '. we can then check that both approaches become very close for 
q <C fcp ; . For higher values of q, the RPA criterion always leads to a reduction of instability compared with the 
semi-classical approach, consistently with the previous discussion of the quantum quenching affecting the dispersion 
relation (Fig. [3]). 
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FIG. 4: (color online) Instability regions associated with fixed values of the wave-number q, at zero temperature, according to 
RPA (thick curves) and semi-classical (thin curves) criterions. Upper part: representation in the density plane. Lower part: 
representation in the plane of Fermi momenta. 



Figure [5] represents the global region of instability against clusterization. For each approach, and for a given 
temperature, it corresponds to the envelope of all the q-spinodals: inside this envelope, nuclear matter can be 
found unstable against finite-size fluctuations for at least one q-value. The right par of Fig. [5] compares the RPA 
instability envelopes for different Skyrme-type forces: the slight differences obtained reflect those appearing between 
the thermodynamic spinodal curves of Fig.O Since the widest q-spinodal contour is obtained for q around 100 MeV/c, 
the corresponding graph of Fig. [4] is representative of the difference between the global instability regions obtained in 
RPA and semi-classical approaches at T = 0. On the left part of Fig. [5] are shown the instability regions obtained at 
T = 5 and 10 MeV for both approaches. The difference between the two envelopes is sensitively enhanced at higher 
temperature, as we approach the limiting temperature of the instability: this means that this limit occurs earlier in 
the RPA due to the quantum quenching. It is then clear that the temperature range concerned by the liquid-gas 
instabilities is not sufficient for a suppression of the quantum effects. 

The temperature dependence of the instability regions is illustrated in Fig. [5] for symmetric matter, considering the 
thermodynamic and finite-size (semi-classical and RPA) spinodals. We can observe the reduction of the thermody- 
namic instability, and the smaller temperature range of the RPA instabilities with respect to the semi-classical ones. 
The temperature and density coordinates of the highest-temperature points appearing on this figure are reported in 
Tab. IIVI for the different Skyrme forces. In the thermodynamic case, the phase-transition region ends in a critical 
point occurring for symmetric matter; its coordinates are (p c ,T c ). In contrast, the finite-size instability region does 
not reach its limiting temperature on the axis of symmetric matter, the charge-exchange symmetry being broken by 
the Coulomb interaction. Instead, the limiting point is given by the coordinates (pu m , {Z/A)y im , Tn m ) (associated with 
the last unstable momentum qn m ). These coordinates are given in Tab.fVl in the semi-classical and RPA cases, for the 
different Skyrme forces. Let us first compare the RPA and semi-classical limiting points. The limiting temperature 
is reduced by several hundreds of keV due to the quantum quenching of the instability in the RPA approach. The 
position of the limiting point is found to be at slightly higher density for the RPA: however, this feature is not very 
pronounced and is not systematic at given Z/A, as can be seen on Fig.|5]for Z/A — 0.5. However, it is noticeable that 
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FIG. 5: (color online) Envelopes of instability against finite-size density fluctuations. Left: comparison between semi-classical 
and RPA criterions, for temperatures T = 0; 5; 10 MeV. Right: comparison between the RPA instability envelopes obtained 
with different Skyrme forces, for T = 0; 10 MeV. 
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FIG. 6: (color online) Total-density limits of the instability region depending on the temperature, along the axis of symmetric 
nuclear matter (Z/A = 0.5). The thermodynamic instability region is represented, as well the finite-size instability regions 
obtained in the semi-classical (SC) and RPA approaches. The symbols indicate the limiting temperature 1f lm in symmetric 
matter: for the thermodynamic curve, it corresponds to the critical point. Left: the three curves obtained with SLy230a. 
Right: comparison of the different Skyrme parameterizations. 



the RPA limiting point occurs at higher asymmetry: this can be attributed to the larger influence of Coulomb effects, 
since, as expected, the corresponding transfered momentum is lower. Let us now compare the critical and limiting 
temperatures. For the forces with similar surface properties (e.g., in the presented results, SLy230a and RATP), 
the difference between the values of Tn m reflects the difference between the thermodynamical critical temperatures 
T c (typically reduced by 2.4 and 2.8 MeV respectively in the semi-classical and RPA cases). A smaller reduction 
of the critical temperature is observed for the force with reduced surface energy (SGII), for which the last unstable 
momentum qi; m is higher. 

We finally explore on Fig. [7] the behavior of the phase separation direction 8pz/8p associated with the unstable 
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pc 

[fm- 3 ] 


s,SC 
Plim 

[fm" 3 ] 


s.RPA 
Plim 

[fm" 3 ] 


T" 1 
-tc 

[MeV] 


T ^sc 
[MeV] 


rriS.RPA 

Aim 

[MeV] 


Sly230a 


0.0538 


0.0461 


0.0458 


14.55 


12.09 


11.63 


SGII 


0.0518 


0.0451 


0.0461 


14.46 


12.33 


11.78 


RATP 


0.0550 


0.0473 


0.0474 


14.72 


12.25 


11.81 



TABLE IV: Limiting density and temperature for symmetric nuclear matter (subscript "s"). The thermodynamic critical 
point is compared with the values obtained from the semi-classical (SC) and RPA approaches. 





Pc 

[fm" 3 ] 


Pfim 

[fm" 3 ] 


RPA 
Plim 

[fm" 3 ] 


(Z/A)fP 


{Z/A)*l A 


[MeV] 


rpSC 

J lim 

[MeV] 


rrnRPA 

-Mini 

[MeV] 


9lim 

[MeV/c] 


RPA 

[MeV/c] 


Sly230a 


0.0538 


0.0467 


0.0474 


0.426 


0.416 


14.55 


12.20 


11.77 


76.7 


70.1 


SGII 


0.0518 


0.0456 


0.0460 


0.447 


0.431 


14.46 


12.39 


11.87 


85.0 


75.1 


RATP 


0.0550 


0.0478 


0.0478 


0.436 


0.427 


14.72 


12.34 


11.92 


76.6 


70.2 



TABLE V: Thermodynamic critical point (total density p c , proton fraction Z/A — 0.5 and temperature T c ) compared with the 
limiting values obtained from semi-classical (SC) and RPA for finite-size instabilities. 

eigen-modes, for discrete values of the momentum q. RPA and hydrodynamics results are compared. The view- 
graphs represented on this figure are separated into two groups: on the left-hand side, the total density is fixed at 
p = 0.05 fm~ 3 and the horizontal axis explores the isospin asymmetry p 3 / p; on the right-hand side, we fix the proton 
fraction Z/A = 0.3 (i.e. pzj p = 0.4), and the horizontal axis explores the total density. In both cases, the instability 
direction is represented for different Skyrme forces (top), different q values (middle), and different temperatures 
(bottom). First, it could be remarked from the left view-graphs that the different curves are mostly situated inside 
the cone defined by the horizontal line and the first diagonal, satisfying the inequalities (|49|) . This confirms the 
predominance of the normal isospin distillation, already noticed about Fig. [3J 

Let us now consider the right panels of Fig. [7] As expected, the difference between RPA and hydrodynamics results 
tends to disappear at the high-density border of the unstable zone, as the Fermi momenta become large enough 
with respect to the transfered momentum q. Indeed, in this case the RPA reduces to the Vlasov approach, which 
identifies with the hydrodynamics on the borders of the instability interval. Inside the spinodal zone, the differences 
are enhanced, especially at T = 0: in the hydrodynamical approach, the instability direction Sp^/Sp strongly decreases 
towards lower densities, while in the RPA approach its density dependence is weak. As a consequence, the RPA values 
of Sps/Sp are larger (up to a factor 2). This accentuates the observation made on Fig. [31 that the distillation effect is 
reduced in the RPA framework. However, both approaches again confirm the isoscalar behavior of the unstable mode. 
We finally notice that the differences are smoothed for the results at finite temperature: this is essentially due to the 
reduction of the instability, limiting the range of directions leading to a negative free-energy curvature, rather than to 
the disappearance of the quantum effects. Indeed, we know from the comparison between RPA and Vlasov curves on 
Fig. [3] that these effects do not affect the separation direction at least for p > 0.05 fin" 3 . Furthermore, as discussed 
previously, quantum effects are still present in all the range of temperature concerned by the liquid-gas instabilities. 

B. Dynamical properties of cluster formation 

In the spinodal decomposition scenario, the early dynamics of cluster formation is induced by the most unstable 
mode, corresponding to the fastest amplification of the fluctuation. For fixed temperature and densities, this is given 
by the maximum of the dispersion relation l/r(q), that we denote (qo, 1/to). Clusters should then be formed according 
to the timescale To, with a size of the order of Xq/2 = 7r/qo, and a composition related to the corresponding direction 
(5p 3 /5p)o. 

We show on Fig. [5] the projection on the density plane of the iso-contours associated with the typical time constant 
To (top) and fluctuation size Ao/2 (bottom). Hydrodynamics and RPA results are shown, respectively on the left 
and right panels. The instability envelope corresponds to To = oo. Inside this envelope, the region of fastest time 
evolution corresponds to the lowest values of To and is found for low-density symmetric matter (p ~ 0.03 fm -3 for both 
approaches). The RPA approach has typical times twice larger than the hydrodynamical ones. The half- wavelength 
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p= 0.05 fm" 3 Z/A = 0.3 
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FIG. 7: (color online) Directions of phase separation for fixed values of the wave number q studied along two different axis of the 
density plane, for the hydrodynamical (thin lines) and RPA (thick lines) approaches. Dots give the limit of instability regions. 
Left: constant density. The horizontal line corresponds to the isoscalar direction, and the first diagonal to the constant-Z/A 
direction. Right: constant proton fraction. Top figures compare different Skyrme forces; middle figures different values of q; 
bottom figures different temperatures. 



Ao/2 is larger in the quantal approach than in the semi-classical one. These results generalize the features observed 
on the dispersion relation l/r(q) illustrated on Fig. [31 i.e. the quantum quenching of instabilities and the higher 
instability in q and 1 /t of the first sound compared to zero sound. 

Let us stress here that, although larger time constants and larger cluster sizes are repeatedly associated in the 
comparison between the different approaches, these two features are actually decoupled. Indeed, tq and Ao are not 
linked by a proportionality relation, since they are determined by the minimum-r relation (dr/dq)o = 0. In fact, the 
observed correlation is due to the particular shapes of the respective dispersion curves l/r(q). It can even be observed 
on the low-density part of the hydrodynamics plots that longer times can be associated with shorter wavelengths, 



18 



CO 



0.06 



0.04 



0.02 



CO 



0.06 



0.04 - 



0.02 L 



Hydro. 



SLy230a, T=0 



RPA 




(X Q /2) levels (fm) 




(V 2 ) levels (fm) 




P (fm 3 ) 

n 



P n Cm -3 ) 



FIG. 8: (color online) Time and space characteristics of the most unstable mode, represented in level-lines inside the instability 
envelope (thick curve). RPA and hydrodynamics results are shown. Top: time constant tq giving the fastest amplification of 
the density fluctuation. Bottom: half- wavelength Ao/2 = 7r/qo of the dominant fluctuation mode. Levels are separated by an 
interval of 0.5 fm. 



stressing the absence of proportionality relation between those two quantities. 

It is interesting to give an estimate of the cluster masses associated with the favored wavelengths represented on 
the bottom part of Fig. [8] Let us first notice that the precise determination of this quantity is beyond the scope of the 
present approach, which describes the initial dynamics of fragment formation; the properties of the fully separated 
fragments in the final state coming from the spinodal decomposition should be obtained following the evolution of 
the density fluctuation beyond the linear regime. However, an order of magnitude can be obtained in the present 
framework assuming that the half- wavelength Ao/2 of the favored mode could give an estimate of the nucleus diameter 
at saturation density. In Tab. IVII are listed some values of the half- wavelength Ao/2 of the favored mode and the 
deduced mass of produced nuclei. Comparing the typical values of the half- wavelength Ao/2 shown on Fig. [5] with the 
masses obtained in Tab. IVII it could be inferred that the spinodal decomposition scenario predicts the formation of 
light nuclei (up to A~ 20) through most of the unstable region. Medium-mass nuclei could be formed in the region 
close to the frontiers of instability (large global density or neutron excess): their formation would thus require much 
more time than that of lighter elements. 

On Fig. [HI we explore the instability direction (Spa/8p)o through the density plane, for RPA and hydrodynamics. 
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FIG. 9: (color online) Phase-separation directions according to the most unstable mode, reprented by levels and tangential lines 
(see text for detailed explanation) inside the instability envelope (thick curve). RPA and hydrodynamics results are presented. 
Levels are separated by intervals of 0.1. 
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16 


25 


37 


53 


72 


96 


125 


159 


199 


244 


masses 


light 


medium 


heavy 



TABLE VI: Relation between the half-wavelength of the favored mode and the mass-number A of a nucleus at normal density, 
A /2 = 2r A 1/3 , with r = 1.2 fm. 



The top view-graphs display level-lines of the quantity {8p^/5p)o- A different kind of representation is employed for 
the bottom view-graphs: they give a direct picture of the instability direction in the plane, by the so-called tangential 
lines. In each point of a tangential line, the tangent to this curve gives the direction of the most unstable mode. This 
representation clearly shows the isoscalar character of the modes obtained in both approaches. 

Let us now comment the level-line panels of Fig. O The (Sp3/6p)o level corresponding to the value is the 
purely isoscalar mode. In the absence of Coulomb interaction, this level strictly corresponds to the axis of symmetric 
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matter |43f , as imposed by the isospin-exchange symmetry. In the present case, we have already noticed that the 
Coulomb interaction disfavors proton-density fluctuations. This results in a global increase of the {5pz/8p)v values, 
thus the levels are shifted towards the proton-rich side. This effect also appears on the left part of Fig. [7] as expected, 
the Coulomb shift is increased for larger values of q (middle graph). Coming back to Fig. let us compare the 
hydrodynamics and RPA upper panels. Two differences can be noticed: i) in RPA the absolute values of (8p3/8p)o 
are larger, and ii) level-lines on the neutron-rich side approximately follow constant-Z/^4 lines while hydrodynamics 
levels are incurved. These panels thus give a global view of the features which have already appeared on Figs. [3] and [71 
namely a reduced distillation in RPA, with a weaker density-dependence of (Sp3/Sp)o at fixed Z/A. 

Figure [10] displays the clusterization properties obtained along selected axis, for different Skyrme forces, at T = 0. 
Time, size and direction characterizing the most unstable mode are addressed, comparing RPA and hydrodynamics 
results. The inverse typical times l/r differ typically by a factor of 2 between the two approaches. The cluster 
sizes are generally increased by about two fermis in the RPA approach. At very low density however, much shorter 
wavelength are favored in the hydrodynamical approach: note that the validity of the Thomas- Fermi approximation is 
under question in this domain. Finally, the curves giving the direction of the most-unstable mode are very similar to 
those obtained on Fig. [Jj for q = 100 MeV/c: this is due to the weak dependence of the instability direction 8pz/8p on 
the momentum q for large enough q values (see Fig. [3]). Let us now compare the different Skyrme forces. Similar results 
are obtained with the three parametrizations. SGII presents faster time evolutions and favors shorter cluster sizes, 
as was already visible on Fig. [3J As for the direction of phase separation, we can notice that the isospin-distillation 
effect is slightly weaker with SLy230a. 

Finally, we explore on Fig. [TTJ the temperature dependence of the clusterization properties. Time, size and direc- 
tion characteristics are shown for T = 0, 5 and 10 MeV, comparing now the three approaches: RPA, Vlasov and 
hydrodynamics. The differences between RPA and Vlasov results give an estimate of the quantum effects. As can 
be seen from the size characteristics, quantum effects are still sensible at T = 10 MeV: the Vlasov curve for qo is 
clearly distinct from the RPA one. The smoothing of the differences of clusterization properties between the three 
approaches with increasing temperature is mainly due to the reduction of the instability. The range of q and Sp3 /5p 
values corresponding to a negative free-energy curvature is narrowed, constraining the different approaches to a closer 
agreement as for the size and direction associated with the most unstable mode. 

The time-evolution part of Fig. [TT1 deserves a special comment. As T increases, all values of 1/tq go towards 
zero, which defines the limit of the instability region. At T = 10 MeV, Vlasov and RPA results for this quantity 
are quite close: however, if we still increase the temperature, the RPA curve will disappear first, then Vlasov and 
hydrodynamics curves will reach zero together. While the 1/to values converge to zero, the corresponding divergence 
of To on the contrary magnifies the differences between the three approaches. As a result, the range of prediction as 
for the time-scale of cluster formation is strongly increased as we approach the limit of the instability region. 

VI. CONCLUSION 

We have presented the clusterization properties of asymmetric nuclear matter in three approaches, comparing the 
quantal RPA results with those obtains in a semi-classical framework, with Vlasov and hydrodynamical calculations. 
In all cases, starting from homogeneous infinite matter, we have studied the consequences of introducing a plane- 
wave density fluctuation of finite momentum q, working in the small-amplitude limit. Liquid-gas instabilities have 
thus been identified, corresponding to an amplification of the fluctuation. The unstable eigen-modes (defined for 
each approach) are characterized by a time constant r, a wavelength A = 27r/q, and a direction dp^/Sp. Assuming 
a scenario of spinodal decomposition, the early dynamics of cluster formation is dominated by the most unstable 
mode (for given temperature and densities). The clusterization properties are then characterized by the quantities 
associated with this mode: a time-scale To, a typical cluster size Xq/2 — 7r/qo, and a cluster composition induced 
by the direction (Sp^/ 8p)$. We have noticed that this study, which assumes fragment formation out of equilibrium, 
is related to the instability properties of homogeneous matter at sub-saturation density; it is complementary to a 
description of equilibrated, clusterized matter. 

Considering infinite nuclear matter is an approximation which permits to treat analytically in a common framework 
the liquid-gas instabilities concerning both compact stars and the nuclear multifragmentation observed in heavy-ion 
collisions. In this framework, the differences between the presented approaches can be attributed to two well-identified 
origins: comparing Vlasov and RPA gives a direct estimate of the quantum effects, while comparing Vlasov and 
hydrodynamics shows the difference between a collisionless and local equilibrium description. 

The results obtained with the different approaches show similar qualitative features, reproduced with the three 
Skyrme-like forces that we have used (SLy230a, SGII, RATP). In particular, the regions of instability against finite- 
size density fluctuations are very close, and the direction of phase separation is always found to be essentially isoscalar, 
favoring the phenomenon of isospin distillation which leads to more symmetric clusters. An order of magnitude of the 




FIG. 10: (color online) Time, size and isospin fractionation characterizing the most unstable mode: dependence on the Skyrmc 
parametrization. RPA (thick line) and hydrodynamics (thin line) results are shown. Left: constant density. Right: constant 
proton fraction. 
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FIG. 11: (color online) Time, size and isospin fractionation characterizing the most unstable mode: temperature evolution. 
RPA (thick line), Vlasov (dashed line) and hydrodynamics (thin line) results are shown. Constant proton fraction: Z/A — 0.5 
for time and size characteristics (four upper lines), Z/A — 0.3 for isospin fractionation. 



cluster masses, estimated as spherical nuclei of diameter Ao/2, indicates that light to medium nuclei (A < 40) should 
be favored through most of the instability region for the various approaches. However, some quantitative differences 
are noticeable. As for the comparison between the Skyrme forces, the main distinction concerns SGII, associated 
with shorter size and time characteristics. Let us now summarize the observations concerning the three different 
approaches. The comparison between RPA and Vlasov shows a quantum quenching leading to larger size and time 
characteristics in the RPA framework. This effect survives for all the temperature range of the instability. The Vlasov 
results are framed between the RPA and hydrodynamics ones, which can be considered as giving two opposite limits 
of predictions. The physical processes should correspond to an intermediate situation, involving different possible 
relaxation rates. This is an important aspect to take into account in realistic calculations aiming at a quantitative 
description of nuclear fragmentation. Indeed the two limits we have obtained in the present work indicate that the 
interval of predictions can be quite large. At T — 0, the inverse time constant 1/tq is typically enhanced by a factor 2 
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between RPA (slower) and hydrodynamics (faster). The cluster size is about 2 fermis smaller for the hydrodynamics, 
with a stronger reduction at low density where the validity of the semi-classical approaches gets more questionable 
(due to the low Fermi momenta). The direction of the unstable modes, for a given value of q, is little sensitive to the 
quantum effects. The main difference is observed between collisionless and local equilibrium approaches: inside the 
instability region, the direction corresponding to a divergence of the response function is slightly less isoscalar than 
the direction of minimal free-energy curvature. 

Approaching the frontiers of the instability region (in densities or temperature), the predictions between all ap- 
proaches get closer as for the size and direction induced by the most unstable mode, due to a reduction of the instability 
interval of q and Sp^/Sp values. The time constant To characterizing the rapidity of cluster formation has a specific 
behavior, since it diverges (by definition) on the borders of the instability region: the range of values obtained with the 
different approaches is then strongly enhanced. This effect is important to notice, since in heavy-ion collisions as well 
as core-collapse supernovae the time-scale of cluster formation should be compared with that of various equilibration 
processes. Conversely, a deeper study of the different time scales involved could indicate the respective validity of 
collisionless and local equilibrium descriptions of matter clusterization. 



APPENDIX A: LINEAR RESPONSE TO AN EXTERNAL PROBE 



We detail the calculation of the response function for the single-fluid case. Let us consider the response to an 
infinitesimal external field of the form: 



V cxt = Ee^-^+^t (t < 0) ; V oxt =0 (t > 0) 



(Al) 



carrying a momentum q and an energy uj (the positive infinitesimal part r\ ensuring the adiabatic switching of the 
external field). At late enough times and small enough external field £, the expectation value of the density fluctuation 
has the same space and time dependence as the external field: 



5p(r,t) = ae 



(A2) 



where a includes the time dependence for a more compact notation. The linear response II (a; , q) is defined by the 
ratio of the density change to the external field strength, for a small external field £: 



I%,q) = - 



(A3) 



A density fluctuation Sp(t) defined such as p(t) = po + Sp(t) induces a variation Sh(t) of the mean-field operator, such 
that h(<) = h + 5h(t). The linearized TDHF equation 



1 8t 



h , <5/3 



reads in momentum space: 



The elements of V cx t are: 



(k|^|k') 



ff(fk' ~ fk) 



Sh + V oxt ,p 



^<k|<m + V ext |k') 



(k|V ext |k') = Se- i ^ + ^ t S kf ^ = e(t)S v+ ^ k 
and the elements of Sp have the form: 

(k\Sp\k') = K (k,q)e- i ^ +il ' )t ^ +q , k = K(k,q,i)^ +q , k 
The density factors a and «(k, q) are related by: 

k 

According to Eqs. (|A5HA7j) . the matrix elements of Sh are of the form: 

<k|5h|k') = ^(k,q)e- i ^+"') t ( 5 k ' +q ,k = ^(k,q,i)^+ q ,k 



(A4) 

(A5) 

(A6) 
(A7) 

(A8) 
(A9) 
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where i/(k, q) is connected to the residual interaction V rcs . Denoting E[/3] the total energy of the system, this residual 
interaction operator is defined by [231 ] : 

5E = E[p + 6p] - E[po\ = V -jr-Spa + ~ Y] a ~ o 5pt 3 $Pki 

= 51 h i^i + 2 XI Y 7i% 5 PiJ 5 Pki 

ij ijkl 

= Tr [Mp] + iTr (1) Tr (2) [v^ 2) ^ (1) ^ (2) ] (A10) 

where we have limited the development to the second order in Spij. From these relations, we also have the expression 
of the mean-field variation in terms of the residual interaction (at first order in 8pij): 

kl F ' 3 kl 

8h = Tr (2) [V^ 2) <5p (2) ] , (A12) 

so that from Eq. (|A9|) we have: 

(k + q|<m|k>=2(k,q,*) = E V k+q,k 3 ;k,k 2 (k 2 |^|k 2 ) 

= E V k+q,k 2 ;k, k2+q ^(k 2 ,q,0 • (A13) 

k 2 

Assuming that the residual interaction clement V k ^ q ka . k k2+q does not depend on (k, k 2 ) - Landau monopolar 
approximation - it can be written: 

V£f q , ka;k)ka+q = V res (q) (A14) 

and we obtain: 

9(k, q, t) = 9(q, t) = V rcs (q) £ R(k, q, t) = W les (q)a(t) = v res (q)5(t) (A15) 

k 

with the notation v ros = fiV ros . The matrix elements of <5h are then: 

<k|<Jh|k') = v res (q)5(t) <5 k , +q . k . (A16) 



Inserting the expressions obtained for the matrix elements of 8p, <5h and V ox t, the linearized TDHF equation (|A5|) 
now reads: 

«(k, q, <)<5 k ' +q ,k = ~ fk \ - . 1 (v res 5(t) + f) i k , +q , k (A17) 

£ s(k, q, t) = na(t) = £ MTT^rr^ i + £ ~) ( A18 ) 

k k ^ - ^kq) + 1??] V / 

where Wkq = Wk+ q — u>k- The second line is obtained by summing over k and k'. Taking the continuous limit of 
Eq. (|A18[) and multiplying by the time factor e _l ' w_1 ' 7 ^ t , we obtain: 



f g(fk - fk+q) ( rcSQ c) 



(27r) 3 - Wkq) + H 

= U Q (uJ,q)(v rcs a + £) . (A19) 

where Ilo is the uncorrelated response function (Lindhard function) . A rearrangement of Eq. (|Al9j) gives the following 
expression for the response function: 

n(^q) = g= "f' q ) ■ (A20) 
£ 1 - n (w,q)v rcs 

This result could also be deduced from the Bethe-Salpeter equation [40| , keeping the transfered- momentum dependence 
of the residual approximation and treating the exchange terms at the Landau approximation: this is the LAFET 
approximation detailed in Ref. [4lj . 
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APPENDIX B: RESIDUAL INTERACTION AND FREE-ENERGY CURVATURE 



In order to establish the link between the residual interaction and the free-energy curvature, we start with the 
thermodynamic case considering a homogeneous system with large particle number TV and volume fl. For simplicity 
we first expose the calculation details in the single-fluid case. Fixing f2, at equilibrium the density matrix p is uniquely 
defined by N and a fluctuation of the particle number SN leads to a fluctuation 5p which can be expressed by the 
Taylor expansion: 



dp 



1 d 2 p 



5/5 = m m+ 2dw 



<5N 2 + o(5N 3 ) 



(Bl) 



and the free-energy variation is: 



<5F b = SE h -TSS 

= ^N+i^N 2 + («5N 3 ) 



(B2) 



where S is the entropy and p the chemical potential. The superscripts "b" (for "bulk") specify that we are in the 
thermodynamic description. i5F b can be separated into an non-interacting term (5F b ' NI and a residual interaction term 
<5F b ' res . The non-interacting part gives the free-energy variation for non-interacting particles of mass m* in a constant 
mean field. It contains the entropy variation, since the entropy S is determined by the density independently of the 
particle interaction. The residual interaction term then reduces to <5F b,rcs = (5E b,rcs , which is the term due to the 
variation of the mean field <51i, so that: 



Referring to the expression (|A10|) . SE h:1 
Taking constant matrix elements for V 



<5F b = <5F b < NI + SE h ' ICS . 
is expressed in terms of the residual-interaction operator as: 



1. 
2 
I 

-Tr (1) Tr (2) 



Tr ( i)Tr (2 ) V^Sp^dp^) 

Ab,res gP(l) ffi(2) 
(1,2) 9N 9N 



<5N 2 +o((5N 3 ) 



b < rcs , such that V]j"Z 



V b < rcs , we obtain (cf Appendix El : 



6E b,re S = Iv b ' rcs <5N 2 + o(6N 3 ) 



(B3) 



(B4) 



(B5) 



so that 



fiSN + - 



NI 



9N 1 



SN 2 + o(<5N 3 



Identifying the second-order term of Eqs. (|B2|1 and (|B6[) gives: 

XI 



dp 



ON 



dp 



^ = l=£) +v 



ON 



/-b,rcs 



(B6) 



(B7) 



The non-interacting part of the chemical-potential derivative is related to the density of states at the Fermi level, 
denoted Nq. Taking the derivative with respect to the density p = AT/O, we have: 



Op 
'dp 



XI 



h 2 dk 2 



2m* dp 

where is the Fermi momentum. Using the notation f2V b ' rcs = v fc 



1 

No 



we can re-express Eq. (|B7I) as: 



dp 
dp~ 



r b,rcs 



(B8) 



(B9) 



which establishes the link between the thermodynamic free-energy curvature and the bulk part of the residual inter- 
action. 
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Let us now consider a finite-size density fluctuation. In order to obtain the expression of the residual interaction 
in terms of the free-energy curvature matrix, we still work in the local equilibrium approach. Following Eq. (|B5p . the 
bulk part of the residual energy variation per unit volume reads (dropping the terms beyond second order in 8p): 

£T^b,rcs 

(BIO) 



9. 



Considering a density fluctuation 5p(r), the average residual energy variation per unit volume is (for the bulk part): 



b , res 



S£ 



With a Fourier development of Sp(r): 



dr 



SE h ' ICS [6p{r)] 




2ft 



dr 



s S P (rf 



A — A* 



(BH) 
(B12) 



we have: 



5£ 



b,rcs 



q 



2 b, res 



(B13) 



Adding a transfered-momentum dependence to the residual interaction, we generalize the previous expression to: 

SS r» = 1 J- |A q | V cs (q) . (B14) 



The q-dependence concerns the gradient and Coulomb terms of the energy variation: 

= I f drH v (r) = U drC* (Vp(r)) 2 = £ |A q | 2 q 2 C v 

J J q 



(B15) 
(B16) 



where ej = 47rg 2 /eo, <?i being the particle electric charge for the considered single-component fluid. We can now 
express v rcs (q): 



5£ rcs = (5£ b < rcs + c5£ v + ^ c = i^|A ql 

q 



r b,rcs 



V„2 



2C v q 



v ros (q) = v b 



2C v q 2 



47re 2 



4ttc 2 



(B17) 
(B18) 



The whole procedure can be easily generalized to the two-component nuclear matter, defining the residual interaction 
operator by the relation (with i,j = n, p): 

d 2 E 



^gres = IV" V" r_ 'V.:k k,'V';:k..k. 

1 



2 £ Tr (l) Tr (2) [V^fj2)^(il)%2) 

Equations (JB9J), (|B10|) . (|B14|) and (|B17|) then become respectively: 



dpi 


b.res , ^'ij 

— v ■ H — 


d P] 


10 N 0>i 


<5E b ' res 



\ A b.rcs c r 


<5£ rcs 


A ■ A* 
;q J;q v 

q i,j 

A ■ A* 
y^ y^ ^WTj'iq 

q i,j 



v b / es + 2C^q 2 + 



(B19) 

(B20) 
(B21) 



(B22) 
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and since 



vJT(q) 



b.rcs 



+ 2C v q2 + 4^ = C f. (q) _ 



N ,i 



(B23) 



we obtain the following relation between the residual interaction and free-energy density curvature matrixes: 

C f = (No)" 1 +v res ; No = ( N J n N ° ) . (B24) 



APPENDIX C: LINK BETWEEN RESIDUAL INTERACTION V ros AND ENERGY VARIATION <5E res 
We first work in a discrete basis |kj), where the residual energy variation is expressed by: 

2<5E rcs = Tr (1) Tr (2) [vft 2) 8p (1) 5p (2) ' 

= E< k i k 2|V res |k 1 k 2 )(k 1 |Mk' 1 )(k 2 |<5pm) (CI) 

In the discrete basis, the one-body density matrix elements are: 

(kilfikj) = nWSij , (C2) 
whose dependence on the particle number N is contained in the occupation number n(kj) = n(kj,p = N/O): 

~dn(ki,p) 



(ki\Sp\kj) = ( ki |^m)5N = i 



dp 



SN . 



We have then for the residual interaction term: 

J®™ 1 V/ k k |yrcs| k h) dn(k u p)dn(k 2 ,p) 



kik 2 



dp 



(C3) 



(C4) 



Let us now take the continuous limit. Denoting (kik 2 |V res |kik 2 ) — ^kTk2'ki k 2 ' ^ c residual interaction term can 
be expressed as: 



<5E ri 



SN 2 



(2 



is y dkid^v^k^k.k. 



5n(ki,p) <9n(k 2 ,p) 



dp 



9p 



(C5) 



We now have to explicit n(k, p). For a fermion gas at thermodynamic equilibrium, at temperature T = 1/(3, with the 
chemical potential p(p) and the individual energies e(k, p), we have: 



n(k, p) = g 



1 + e 



-/3X(k,p) 



X(k,p) = /i(p)-e(k,p) 



(C6) 



with e(k, p) the individual particle level of momentum k and g the level degeneracy. The residual interaction term 
now reads: 



5W 



1 f AVAV V- dX{k u p) dn(X{k u p)) 8X{k 2 , p) dn(X(k 2 ,p)) 

dkidk 2 V kl k2 . kl k2 



SN 2 (2tt) 6 J ~" L ~ " a p M 

For non-relativistic nucleons, with a Skyrme interaction, we have: 

fc 2 



dp 



<9X 



X(k,p)=X(/^,p) = p(p) 



U(p) 



2m* (p) 



Hp) 



k 



2m* (p) ' 



At zero temperature, the density derivative of the occupation number is a Dirac function: 

n(fc 2 ,p) = gQ{X{k\p)) 



dn(k 2 ,p) 
dp 



= 9- 



d X(k 2 ,p) 

dp 



5{X{k\p)) 



(C7) 

(C8) 

(C9) 
(CIO) 
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Since X = for k 2 = fc| , the integral becomes 



gk F m* dX(k 2 ,p) 
(2tt) 2 dp 



J dcos#idcos 2 Vp s (cos6' 1 , cos 9 2 ) 



(Cll) 



Neglecting the angle dependence of the residual interaction, we have V F cs (cos 9\, cos 62) — Vp cs . Using the relation 

dX(p,k 2 ) _ 2ir 2 



we finally obtain that 



dp 



SE 1 



m*gfcj 



SN 2 



(C12) 



(C13) 



In order to generalize this result at finite temperature, let us now assume that the residual interaction elements 



"\rres 

v ki,k 2 ;ki,k 2 



are equal to a constant V res . Starting again from Eq. (|C5[) we simply have 



2^: = v- 1 



SN 2 



(2tt) 6 



dk ld k / n(kl ' p)an(k2 ' p) =v- 



dp 



dp 







dk 



1 2 



dp J (2tt) 3 



n(k, p) 



V r 



(C14) 



whatever the form of n(k, p), since p = J / 27! \3 n(k, p). 
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